Coulomb tunneling for fusion reactions in dense matter: 
Path integral Monte Carlo versus mean field 
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We compare Path Integral Monte Carlo calculations by Militzer and Pollock (Phys. Rev. B 71, 
134303, 2005) of Coulomb tunneling in nuclear reactions in dense matter to semiclassical calculations 
assuming WKB Coulomb barrier penetration through the radial mean-field potential. We find a very 
good agreement of two approaches at temperatures higher than ~ ^ of the ion plasma temperature. 
We obtain a simple parameterization of the mean field potential and of the respective reaction rates. 
We analyze Gamow-peak energies of reacting ions in various reaction regimes and discuss theoretical 
uncertainties of nuclear reaction rates taking carbon burning in dense stellar matter as an example. 

PACS numbers: 



I. INTRODUCTION 

Nuclear fusion reactions in dense stellar matter affect 
the evolution of ordinary stars and compact stars such 
as white dwarfs and neutron stars. Hydrogen and he- 
lium burning, and later the burning of carbon and heav- 
ier elements |l| drives an ordinary star through the main 
sequence and giant /red-giant branch towards its final mo- 
ments as a normal star. Explosive burning of carbon and 
other elements in the cores of massive white dwarfs trig- 
gers type la supernova explosions (see, e.g., [2] and refer- 
ences therein). Thermonuclear burning of accreted mat- 
ter in surface layers of neutron stars, which enter compact 
binaries, produces type I X-ray bursts ^3]. Deeper burn- 
ing of carbon in accreting neutron stars is likely respon- 
sible for superbursts observed from some X-ray bursters 
(e.g., Refs. Even deeper burning of accreted mat- 

ter in pycnonuclear reactions in the crust of transiently 
accreting neutron stars can power thermal radiation ob- 
served from neutron stars in soft X-ray transients in qui- 
escent states (see, e.g., Refs. [1, @, AH in all, nuclear 
fusion is important in all stars at all evolutionary stages. 

It is well known that nuclear reaction rates in dense 
matter are determined by astrophysical 5-factors, which 
characterize nuclear interaction of fusing atomic nuclei, 
and by Coulomb barrier penetration preceding the nu- 
clear interaction. We will mostly focus on the Coulomb 
barrier penetration problem. Fusion reactions in ordi- 
nary stars proceed in the so called classical thermonuclear 
regime in which ions (atomic nuclei) constitute nearly 
ideal Boltzmann gas. In this case the Coulomb barrier 
between reacting nuclei is almost unaffected by plasma 
screening effects produced by neighboring plasma par- 
ticles. The Coulomb barrier penetrability is then well 
defined. 

However, in dense matter of white dwarf cores and 



neutron star envelopes the ions form a strongly non-ideal 
Coulomb plasma, where the plasma screening effects are 
very strong. Plasma screening greatly influences the bar- 
rier penetrability and the reaction rates. Depending on 
the density and temperature of the matter, nuclear burn- 
ing can proceed in four other regimes [8|. They are 
the thermonuclear regime with strong plasma screening, 
the intermediate thermo-pycnonuclear regime, the ther- 
mally enhanced pycnonuclear regime, and the pycnonu- 
clear zero-temperature regime. The reaction regimes will 
be briefly discussed in Sec. [TTl In these four regimes the 
calculation of the Coulomb barrier penetration is a com- 
plicated problem. There have been many attempts to 
solve this problem using several techniques but the exact 
solution is still a subject of debates. Various techniques 
and approaches will also be outlined in Sec. [Ill 

In Sec. lIIII we analyze recent Path Integral Monte Carlo 
(PIMC) calculations of fusion reaction rates by Militzer 
and Pollock [9] . We compare these results with those ob- 
tained within a much simpler formalism of semi-classical 
Coulomb tunneling in a mean-field potential. To this 
aim, in Sec. IIVI we analyze and parameterize the mean 
field potential in a strongly non-ideal classical ion plasma 
and also calculate and parameterize mean-field reaction 
rates. Section|V]is devoted to comparison of calculations 
by different authors. We conclude in Sec. IVII 



II. NUCLEAR REACTION REGIMES 

A. Physical parameters 

We consider a fusion reaction between identical nuclei 
{A, Z) -\- (A, Z) in a one-component plasma (OCP) of 
atomic nuclei (ions) in dense matter. Here, A and Z are 
the mass and charge numbers of the nuclei, respectively. 
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FIG. 1: (Color online) Temperature-density diagram for car- 
bon matter; T; (T — 1) is the temperature below which ions 
form strongly coupled liquid; Tp « 0.513) is the ion plasma 
temperature; Tm is the solidification temperature in a classi- 
cal ion liquid; Tq is the temperature below which the carbon 
burning rate is temperature independent; I-V label domains 
of different nuclear burning regimes (Table nil. The shaded 
region is most important for applications of carbon burning; 
it is restricted by the lines of constant burning times (Is and 
10^^° yr, for the upper and lower lines, respectively, from Ref. 
[lo|). Filled dots show some T — p points for which PIMC cal- 
culations of reaction rates have been performed if applied 
to carbon burning. 



The results will be general but they will be illustrated 
taking the ^^C+^^C reaction as an example that is most 
important for astrophysical implications (Sec. The 
temperature-density (T — p) diagram for carbon matter 
is shown in Fig.[TJ The filled dots show some T — p points 
for which PIMC calculations [oj have been performed if 
applied to carbon burning. Notice that the majority of 
the data points correspond to much higher T and 
p, where carbon would be immediately transformed into 
other elements either through beta captures or through 
intense nuclear reactions. All the PIMC data @ are ana- 
lyzed in Secs. lIIIflVl The shaded region in Fig.[l]is briefly 
described in Sec. Ill Gl 

Under the conditions displayed in Fig.[Tl carbon is fully 
ionized and immersed in a nearly uniform electron back- 
ground; the electrons are mostly strongly degenerate. 

Coulomb coupling of ions is determined by the familiar 
parameter 



F = Z^eyikBaT), 



(1) 



where a = [3/(47rni)]^/'^ is the ion-sphere radius, Ui is the 
ion number density, and ks the Boltzmann constant. At 



TABLE I: Nuclear reaction regimes in dense matter 



Line 



Regime 



Domain 



I Thermonuclear with weak screening 

II Thermonuclear with strong screening 

III Thermo-pycnonuclear 

IV Thermally enhanced pycnonuclear 

V T = pycnonuclear 



Tp<T<Ti 
0.5Tp <T<Tp 
Tq<T < 0.5Tp 



temperatures T :s> Ti = Z'^e^ /{kBo) (i.e., at F < 1) the 
ions constitute an almost ideal Boltzmann gas. At T < T; 
(r ^ 1) they smoothly (without any phase transition) 
transform into a strongly coupled liquid. The ions solid- 
ify at T = Tm which is much lower than Tj . The freezing 
into the classical body-centered cubic (bcc) crystal oc- 
curs at Tm ~ 175, that is at Tm — Ti/Tm- The difference 
of free energies of various Coulomb structures at low tem- 
peratures is very small, and the actual microstructure at 
these temperatures is rather uncertain (as discussed, e.g., 
in EIHI). 

The importance of quantum effects in motion of plasma 
ions can be characterized by the ion plasma temperature 
Tp determined by the ion plasma frequency ujp, 



hwp/kB, ujp = \/A7rZ^e'^ni/mi, 



(2) 



where is the ion mass. At T >Tp quantum effects are 
relatively weak and quantization of ion motion is mainly 
unimportant. At lower T quantum effects become most 
essential. At p > 4 x 10^^ g cm^"^ dense matter contains 
free neutrons dripped off atomic nuclei. 

OCP of ions can also be characterized by the parame- 
ters 



as 



4F2 



1/3 



37r2T2 



1/3 



(3) 

where qb = / {Z'^e^rrii) is the ion Bohr radius. A state 
of ions is determined by two parameters, for instance, 
by r and T/Tp (or C), while other parameters can be 
expressed through these ones. 

We will also need a temperature (Fig. [1]) 



Tq^Q.bTp/\n{Ti/Tp), 



(4) 



which is the upper temperature of temperature- 
independent pycnonuclear burning (Sec. Ill Pjl . 

The five nuclear reaction regimes Q are summarized in 
Table IT] Temperature-density domains for carbon burn- 
ing in these regimes are seen from Fig. [TJ We outline 
these regimes for fusion reactions of identical nuclei in 
an OCP. 
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B. Thermonuclear burning with weak plasma 
screening 

This regime (regime I in Table |l| is realized at T ^ 
Ti. In this case Coulomb coupling of ions is weak and 
Coulomb tunneling is only slightly affected by plasma 
screening effects. The main contribution to reaction rate 
comes from a small amount of ions with energies E near 
the Gamow-peak energy i?pk (much higher than kBT). 

Neglecting the plasma screening effects, one obtains 
the familiar thermonuclear reaction rate 



th 



'2-Epk ^(-Epk) 
3/i ksT 



exp(-r), 



(5) 



where = ksTr/S, /i = TOi/2 is the reduced mass, 
S{E) is the astrophysical factor (assumed to be a slowly 
varying function of E), and 



1/3 



(6) 



The plasma screening enhances the reaction rates. In 
the thermonuclear regimes with weak and strong screen- 
ing (at T > Tp, and actually even at somewhat lower T) 
the "screened" rate can be conveniently written as 



pscr 
^th 



th -'^^scr ; 



exp(/i). 



(7) 



where i?th is given by Eq. ([5]) and -Fscr is the enhance- 
ment factor expressed through a function h. In the weak 
screening Debye-Hiickel limit (F ^ 1, T ^ T;), one ob- 
tains [III /i = V^T^/"^ < 1, so that Fscr « l-t- /i is only 
slightly higher than 1. 



C. Thermonuclear burning with strong plasma 
screening 

This regime (regime II in Table H]) occurs at Tp < T < 
Ti . The majority of ions are strongly coupled by Coulomb 
forces in their potential wells; quantum effects in their 
motion are weak. The main contribution to the reac- 
tion rate comes from a small amount of highly energetic 
ions which are nearly free and have Gamow-peak energies 
(modified by the screening effects) . The plasma screening 
strongly enhances the reaction rate. 

The screening is often modeled assuming that the re- 
acting nuclei move in a potential 



U{r) = Z-'eZ/r - H{r), 



(8) 



where H{r) is a static and spherically symmetric mean- 
field plasma potential. This approach neglects fluctua- 
tions of plasma screening microfields during an individual 
tunneling event. 

Generally, the function h can be split into two parts. 



The leading term ho = H{0)/kBT is calculated assum- 
ing a constant mean-field plasma potential H{r) = i?(0) 
during the quantum tunneling, while hi is a correction 
owing to a weak variation of H(r) along the tunneling 
path and owing to possible deviations from the mean-field 
approximation (concerned with fluctuations). Estimates 
show (e.g., Ref. [l^) that typical tunneling lengths of 
the reacting ions in the thermonuclear regime (T > Tp) 
are smaller than the ion sphere radius a, and typical 
tunneling times are shorter than the plasma oscillation 
time scales ~ c^p^- This justifies the approach of almost 
constant and static plasma potential during a tunneling 
event as a first approximation. 

The main screening quantity Hq is determined by H{0). 
For a classical ion system, H{0) can be calculated as 
H{0) = AT, where AT is a difference of the Coulomb 
free energy for a given system of nuclei and for a system 
with two nuclei merged into one compound nucleus (e.g., 
DeWitt et al. [l3|)- In this case the leading enhancement 
factor exp(ft,o) — exp(AJ^/fcBr) depends on the one ar- 
gument F. 

The mean- field potential H{r) for a classical strongly 
coupled OCP of ions (liquid or solid) can be determined 
from classical Monte Carlo (MC) sampling (e.g., De- 
Witt et al. |l4j). We will analyze the latest results in 
Sec. IIV Al MC sampling gives the static classical ra- 
dial pair distribution function of ions g{r) which equals 
g{r) = exTp[—U{r)/kBT]. In this way one obtains accu- 
rate values of g{r) and H{r) ~ Z'^e? jr -f kBT\ng{r) at 
not too small r (typically, at r > a), because MC statis- 
tics of close ion separations r < a is poor due to strong 
Coulomb repulsion of ions at small distances. The po- 
tential H(r) at small r, required for calculating the tun- 
neling probability, can be obtained by extrapolating MC 
values of H{r) to r ^ 0. The extrapolation is a delicate 
procedure (as shown, e.g., by Rosenfeld [Tsjl. 

Assuming a linear mixing rule in a multi-component 
strongly coupled ion plasma, Jancovici .l&j obtained 



/io(r) = 2/o(r)-/o(2^/3r). 



(10) 



ho + hi, Fscr = exp(/io) exp(ft,i) 



(9) 



where /o(r) is a Coulomb free energy per one ion in an 
OCP (in units of fcsT). Using MC data available by that 
time (1977) he got ho{r) given in line (a) of TableHIl The 
same expression was used by Itoh et al. [T^ . 

Recent MC calculations for a classical Coulomb liquid 
at F > 1 give highly accurate values of /o(F) (accurately 
approximated by analytical functions, e.g., Ref. ^3]) and 
confirm the validity of the linear mixing rule (e.g., Ref. 
^]\). The function ho{T) has been calculated from Eq. 
PU)) in many papers (e.g., [H, [H, [l^ Ull), and the re- 
sults are in good agreement. In line (b) of Table |TT] we 
present an analytical approximation of /io(r)j which fol- 
lows from the recent MC results of DeWitt and Slattery 
[2^ for a Coulomb liquid. It seems to be the best avail- 
able evaluation of holT), which is in very good agree- 
ment with the expression in line (a). In the indicated 
F-range, it is accurately approximated by a linear func- 
tion /10(F) « 0.9(2^/3 - 2) F = 1.0573 F suggested by 
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TABLE II: Function ho{T) as calculated by different authors. 

Line Ref. feo(r) T 

(a) Eq. (17) in [16] 1.0531 T + 2.2931 T^/* - 0.5551 In F - 2.35 1 < T < 155 

(b) Eq. (20) in [10] 1.0563 T + 1.0208 r"'^^^^ - 0.2748 In T - 1.0843 1 < T < 170 

(c) Eq. (6) in [17] 1.148 L - 0.00944 T In T - 0.000168 r(ln T)^ 5 < T < 180 

(d) Eq. (19) in [18] 1.132 L - 0.0094 Tin T 1 < T < 170 



Salpeter using a simple ion-sphere model. At F ^ 1 
the plasma screening enhancement is huge. For instance, 
exp(/io) ~ 10^'' for F ~ 170. 

Some authors calculated H{0) and the related enhance- 
ment factor exp(/io) by extrapolating MC H{r) to r ^ 
(as mentioned above). In particular, Ogata et al. [l7l[23| 
used that formalism to analyze the enhancement of nu- 
clear reactions in one-component and two-component 
strongly coupled ion liquids. Their /io(r) is given in line 
(c) of Table [TTl These calculations are less accurate than 
those based on Eq. pO]) because of theproblems of ex- 
trapolation of H{r) to r ^ in Refs. [13, [ll] (see Refs. 
[ia mill for details). 

The function /io(r) was also calculated by Ogata [l^ 
using the PIMC method. His result [line (d) in Table [H] 
is in better agreement with the most accurate result [line 
(b)] as discussed in Ref. [l0| . 

In addition to exp(/io)j the enhancement factor i^scr 
in Eq. ^ contains a factor exp(/ii) which depends on 
two arguments (e.g., F and Q). The basic term in hi 
in the thermonuclear regime with strong screening was 
obtained by Jancovici [iq : it is presented in line (A) of 
Table IIIIl where C, is given by Eq. ^ . Introducing t 
from Eq. (O we have C = 3F/t « n/a ~ (Tp/T f^^, n 
being the Coulomb tunneling length in the thermonuclear 
regime (T > Tp). Clearly, ( can be regarded as a small 
parameter in that regime. This basic term can be eas- 
ily obtained in the mean-field approximation (Sec. IIVB[) 
with the lowest-order expansion terms of H(r) over r, 
H{r) = H{0) - \ [rjaf Z'^e^/a. Treating the correc- 
tion as small and using the semi-classical approximation 
for the tunneling probability, one immediately comes to 
Eq. (A). 

Equation (A) in Table IIIIl can be treated as the well 
defined lowest-order term in the expansion of /ii(F, <^) in 
powers of C- There were several attempts to improve Eq. 
(A) by adding new terms obtained either on theoretical 
grounds or by fitting numerical results. These new terms 
are model dependent and debatable. It is thought that 
adding these terms allows one to extend the results to 
lower temperatures, somewhat beyond the lowest bound- 
ary T ~ Tp for the thermonuclear regime (let us remark 
that T ^ Tp corresponds to C = 0.513, and T = Tp/2 
corresponds to ( = 0.815). We will discuss the validity 
of such extensions in Sec. |Vl 

Alastuey and Jancovici Q proposed semi-analytic 
corrections to (A). Their result is given in line (B) of 



Table UTIl Ogata et al. ^ HI] calculated the Coulomb 
tunneling probability and /ii(F, C) using the mean-field 
potential and solving numerically an effective radial 
Schrodinger equation. A fit to their calculations is given 
in line (C). These results were used by Kitamura [1^ 
for constructing an analytic expression for nuclear re- 
action rates in all burning regimes. Ogata [18] calcu- 
lated ft,i(F, C) using PIMC. His fit is presented in line 
(D). Itoh et al. 13] determined the enhancement factors 
Fscr calculating the WKB Coulomb barrier penetrability 
in a mean-field potential. Their results are equivalent to 
hi{r, C) = 1.25 F - r /(C) - ho{r), where /io(F) is given 
by Eq. (a) of Table [III and f{() is given by their lengthy 
fit expression (4.4) [with our ( denoted by /3 and ( < 5.4]. 
Let us remark that their mean- field potential H{r) is sim- 
plified. Its H{0) value is correct but small-r behavior is 
approximate (a continuous function with a break). 

Recent PIMC calculations by Militzer and Pollock Q 
will be analyzed in Sec. El 

D. Zero-temperature pycnonuclear burning 

Zero-temperature pycnonuclear regime (regime V in 
Table H]) takes place at low temperatures, T <Tq, partic- 
ularly, at T = 0; Tq is defined by Eq. ([4]) and plotted in 
Fig-Ill lu this case, one can safely assume that all plasma 
ions occupy ground states in their potential wells. Quan- 
tum tunneling and nuclear fusion occur mainly between 
close neighbors owing to zero-point ion vibrations; ther- 
mal effects in ion motion are unimportant and the reac- 
tion rate is temperature-independent. The reaction rate 
increases with growing p because zero-point vibrations 
become more efficient. 

Coulomb tunneling probability in this regime has been 
calculated in various approximations (see Refs. [l^, |Tl| 
for a recent analysis of the results). The predicted re- 
action rates have similar density dependence but differ 
within several orders of magnitude. 

E. Thermally enhanced pycnonuclear burning 

Thermally enhanced pycnonuclear regime (regime IV 
from Table HI occurs at higher temperatures, Tq < T < 
0.5 Tp. These temperatures are still so low that the ma- 
jority of ions occupy their ground states. However, some 
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TABLE III: Function hi{r,() as calculated by different authors. 



Line Ref. fei(r,C) T C T/Tp 



(A) 


Eq. (35) 


in 


[16J 


-(5/32) rc^ 


1 < r < 155 


C< 1 


T/Tp > 0.37 


(B) 


Eq. (28) 


in 


[24] 


-(5/32) r + o.oi4rc^ + 0.0128 rc'' 


1 < r < 155 


C < 1-6 


T/Tp > 0.18 


(C) 


Eq. (36) 


in 


[17] 


-(5/32) [1 + (1.1858 - 0.2472 logL) C - 0.07009C^] 


5 < r < 180 


C<2 


T/Tp > 0.13 


(D) 


Eq. (19) 


in 


m 


-(5/32) r (2(1 - 0.0348C - 0.1388(2 + 0.0222C^) 


1 < r < 170 


C<2 


T/Tp > 0.13 



of them populate higher bound energy levels in the po- 
tential wells. These ions give the major contribution to 
the reaction rate because it is much easier for them to 
penetrate through the Coulomb barrier. 

This burning regime was studied in Refs. [1, [11] ■ As 
discussed in Ref. Yu\, the results differ within several 
orders of magnitude and have to be improved. 

When the density grows up, the upper and lower 
boundary temperatures for this regime, O.bTp and Tq, 
become closer (Fig. [IJ and finally merge. Therefore, this 
reaction regime disappears at sufficiently high densities 
(when spacings ~ fuvp between quantum energy levels in 
potential wells become too large). As a rule, such a den- 
sity is too high to be of practical importance (it would 
be ~ 10^^ g cm~^ for carbon burning, in which case no 
carbon can survive in dense matter). 

F. Thermo-pycnonuclear burning 

This regime (regime III from Table H]) takes place at 
Tp/2 < r < Tp. It is intermediate between thermonu- 
clear regimes and pycnonuclear ones. When the temper- 
ature increases from ~ Tp/2 to ^ Tp, those ions, which 
give the most important contribution into the reaction 
rate, become unbound and move to continuum states 
(from closest-neighbor collisions in the pycnonuclear case 
to collisions of freely moving particles in the thermonu- 
clear case). 

G. All regimes and problems 

Several works [13, [U, [IB] have suggested analytic fits 
for the nuclear reaction rates valid in all five burning 
reg imes. In particular, the fits constructed in Refs. |10l . 

take into account current theoretical uncertainties of 
the reaction rates and give the optimal, maximum and 
minimum theoretical rates. 

In Fig. [1] we present the lines which divide the T — p 
diagram into five regions appropriate to five regimes. In 
addition, we plot two lines [lO| along which the character- 
istic carbon burning times Tburn = n-i/R are equal to 1 s 
and 10"'^'' years (the upper and lower lines, respectively). 
The astrophysical factor S{E) for the carbon burning is 
taken from Ref. (We neglect a possible hindrance 



of the carbon reaction [27| at energies much lower than 
the Coulomb barrier energy.) The lines are almost hor- 
izontal at lower p and higher T, where carbon burns in 
thermonuclear regimes, and they are almost vertical at 
higher p and lower T, where carbon burns in pycnonu- 
clear regimes. Above and to the right of the Tburn = 1 s 
line the burning is very fast, and there is no carbon in 
dense matter. Below and to the left of the Tburn = 10^" yr 
line the burning is extremely slow (practically absent). 
Therefore, the density-temperature domain of practical 
interest for carbon burning is located between these two 
lines. 

Although the main features of the Coulomb tunneling 
in dense matter seem clear, some important tunneling 
problems in OCP are still unsolved: 

(i) In the thermonuclear regime, where the plasma 

screening is conveniently described by the enhance- 
ment factor Fscr [Eqs. ([7]) and the function 
ho{T) is well defined [line (b) of Table HI]. But 
what is an exact form of the function hi{T, () (Ta- 
ble [ml? 

(ii) Down to which temperatures this description can be 

extended? 

(iii) What is an exact expression for the reaction rate 
at lower temperatures when the burning becomes 
pycnonuclear? 

Our aim will be to analyse problems (i) and (ii) and dis- 
cuss (iii) . The problems in multicomponent ion mixtures 
are more complicated [TT| . 

III. MILITZER-POLLOCK PIMC 
CALCULATIONS 

New PIMC calculations of Militzer and Pollock ^ were 
performed for a wide range of plasma parameters which 
cover all regimes of nuclear burning. In Fig. [1] filled dots 
show some densities and temperatures of stellar matter 
which correspond to four Militzer-PoUock points (the val- 
ues of r and Vs in their Table 1) if applied for carbon 
burning. The majority of their points are not shown; 
they would refer to much higher T and p than those dis- 
played in Fig. [T] (far from the T — p region where carbon 
can exist in dense stellar matter). 
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In principle, the PIMC is the best method to calculate 
Coulomb tunneling in nuclear reactions. It can take into 
account all the effects of dense environment on Coulomb 
tunneling, including fluctuative nature of plasma poten- 
tial; dynamical response of plasma ions to the motion of 
reacting nuclei in the course of quantum tunneling; fi- 
nite width of trajectories of the tunneling nuclei (that is 
beyond the WKB approach). Unfortunately, highly accu- 
rate PIMC calculations require huge computer resources 
(long PIMC runs with many plasma ions involved) while 
all PIMC calculations performed so far are naturally lim- 
ited, at least by not too many plasma ions. 

Notice that the PIMC simulations @ neglected the 
effects of quantum statistics of ions. As shown by several 
authors (e.g., Ref. [3), this is a good approximation 
for the conditions of practical interest. Usually, nuclear 
burning occurs when p is still insufficiently high for the 
quantum statistics effects to be pronounced. 

PIMC simulations of Militzer and Pollock [3] included 
54 plasma ions. The authors calculated the contact 
probabilities g(0), which are the values of the quantum- 
mechanical radial pair distribution function g(r) at r — > 
0. For OCP of ions, g{0) is related to the reaction rate 
through [H] 



i?=^^5(i^pk)5(0), 
TT n 



(11) 



where as — H'^ /{rmZ^e^). 

Militzer and Pollock present their g{0) = gMp{0) for 
36 values of plasma parameters F and rj (equivalently, 
for 36 values of p and T, four of which are shown in 
Fig. [1]). Specifically, they considered 10 values of F: 
F = 0.5 (five ry-points), F =1 (five points), F —2 (five 
points), F =5 (five points), F =10 (five points), F =40 
(five points), F =100 (three points), F =200 (one point), 
F =400 (one point), and F =600 (one point). All their 
results are shown in Fig. [21 where we compare them with 
our calculations in the mean-field WKB approximation 
[5(0) = gMF(O); see Sees. |IV] and |V] below]. Note five 
typos in Table I of Ref. [9j]; in the values of — ln[g(0)] in 
column 5 for rj = 0.25, F = 100; r] = 0.5, F = 100; and 
1] = 1, F=200, 400 and 600, one should remove zero after 
dot; this removal restores correct errobars of — ln[g(0)] 
(otherwise the errorbars are ten times smaller than their 
actual values). 

We have divided the Militzer-PoUock data (somewhat 
arbitrarily) into three groups (i), (ii) and (iii) (listed in 
Table |IV|. 

The first group (i) consists of five points with F = 0.5 
(marked by squares in Fig. [2]). These points correspond 
to a moderately strong Coulomb coupling (T = 2T;, 
thermonuclear burning intermediate between weak and 
strong plasma screening regimes). We do not analyze 
them in detail because in this case the effects of plasma 
screening on Coulomb tunneling are weak and exact 
screening enhancement factors have not been calculated 
so far by other theoretical methods. 
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FIG. 2: Contact probabilities versus F for five values of r;=0.1, 
0.25, 0.5, 1, and 2. Left vertical axis: Filled symbols show 
all Militzer-PoUock gMp{0) data; squares, dots and triangles 
mark the data of the three groups (i)-(iii) (Table lrV)l . Solid 
lines are the contact probabilities (;mf(0) calculated in the 
mean-field WKB approximation. Right vertical axis: Open 
symbols with errorbars display the ratio gMp{0)/ gMF{0) of 
the PIMC to mean-field WKB results [dotted lines refer to 
5mp(0) = 0mf(O) to guide the eye]. Dashed lines show the 
ratio SmV(0)/5mf(0) of our mean-field WKB fitted fSec. UVCl) 
and calculated values. See text for details. 
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TABLE IV: Three groups of Militzer- Pollock data points [|]. 



Line 



Group 



Domain 



(i) Moderate plasma screening F — 0.5 

(ii) T dependent rate, strong screening F > 1, C < 3 

(iii) T independent pycnonuclear data (see text) 



The second group (ii) includes 26 data points (dots in 
Fig. ^ with 1 < r < 100 and C < 3 (T > 0.07irp). 
Three such points belong to the thermonuclear regime 
with strong screening (T > Tp), seven points corre- 
spond to the intermediate thcrmo-pycnonuclear regime 
(0.5Tp < T < Tp), while other 16 points refer to pycnonu- 
clear burning at not too low T (0.071Tp < T < 0.5Tp). 

The third group (iii) includes six points (triangles in 
Fig. HI). Three of them have very large F = 200, 400 
and 600 (with C = 4.33, 5.45, and 6.24, respectively). 
The other three have lower F = 40, 40, and 100 but 
large C (2.53, 3.18, and 4.32), i.e., smaU T < 0.0912rp 
at which the contact probability is expected to become 
temperature independent. The point with F = 40 and 
C = 2.53 belongs also to the second group. 
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IV. MEAN FIELD WKB APPROXIMATION 

In this section we calculate the nuclear reaction rate in 
a mean-field WKB approximation, which is much simpler 
than the PIMC. 



FIG. 3: (Color online) Normalized mean field MC potential 
u{x) in an OCP of ions for 12 values of F from 1 to 800. The 
data for F < 200 are obtained for ion liquid, while the data 
for F = 400 and 800 correspond to bcc ion crystal. Filled and 
open dots show data points included into fitting and excluded 
from it, respectively. Lines show the fit; A'' is the number of 
data points included into the fitting. See text for details. 



A. Mean-field potential 

We start with the discussion of the mean field potential 
H{r) in a strongly coupled OCP. We have taken the re- 
sults of extensive Monte Carlo (MC) calculations of H{r) 
for a classical OCP of ions in a liquid phase for ten val- 
ues of F from 1 to 200 within some intervals of x — r/a 
(Fig. [3]) . Although a thermodynamically stable phase of 
OCP at F > 175 is a crystalline bcc lattice, the differ- 
ence of free energies of the lattice and liquid is small. As 
a result, the OCP can stay liquid at temperatures much 
below the melting point, and one can simulate super- 
cooled liquid in MC runs. In order to check our results 
we have also taken MC calculations of H{r) for a clas- 
sical bcc Coulomb crystal at F = 400 and 800 (two last 
panels in Fig. [3|). 

It is useful to introduce a dimensionless mean-field po- 
tential u{r), 

H{r)/kBT ^Tu{x), (12) 

which is plotted in Fig. [3] 

The dimensionless potential u{x) depends on two vari- 
ables, X and F. For further use, we have fitted the calcu- 



lated values of u{x) by an analytical expression. The fit- 
ting included selected MC points (denoted by filled dots) 
for an OCP liquid at 1 < F < 200. The number N of 
these data points at any F is presented in Fig. [3] (for a 
better visualization we have not plotted some points). 
Other MC points, not included into the fitting, are de- 
noted by open dots. In particular, we did not fit the data 
for the ion crystal (F = 400 and 800). We have also ex- 
cluded from the fitting several data points at low x for 
F < 5 and F = 160. These points were calculated with 
large errors because of poor MC statistics at low sepa- 
rations x. We do not show their errorbars to simplify 
the figure. The errorbars for all other points are small 
and would be invisible. The data points at sufficiently 
large x for F < 200, excluded from the fitting, have been 
used for checking the quality of our fit. Similar poten- 
tials for OCP liquid have been computed, e.g., in Rcf. 
[ii]. The potentials suggested in Ref. [13] for 4 < F < 90 
and a; < 1.6 are also very accurate. 

To ensure a high fit accuracy at small where MC 
statistics is insufficiently good, we have taken into ac- 
count that at small x the function u{x) can be expanded 
in powers of x. In the limit of strong Coulomb coupling 



8 



(r ^ 1) only even powers of x should survive [3l| and 
the expansion should have the form 



u{x) = af) + a2 x'^ + 0:4 a;^ + x^ + 



(13) 



where ao = /io(r)/r, a2 = —j UM, and 04, ag, . . .can 
depend on T. Thus, ao is a slowly varying function of 
r determined by the function ho{T) (discussed in Sec. 
IIIC|) . One can use any accurate representation of ho(r), 
for instance, Eqs. (a) or (b) from Table |TT1 However, 
the fit formula (fH)) for u{x) presented below is especially 
accurate if hoiT) is given by our own fit expression p9p 
(see Sec. IIV C|) . For a not too strong Coulomb coupling 
(r < 5), odd powers of x can become pronounced in the 
expansion (|13p . and a2 can be modified. 

We have fitted our MC data points by the analytic 
expression 



u{x) = ao 
where 



1 - C4 a; - 2 (Ci/ao) x^ + C3 a;" + C2 a;^ 



1 + C2 ag a; 



2 ,,10 



Ci = 0.25 - 0.267r- 



0.084- 0.144 T" 



C2 = 0.05, 
d = 0.434 T" 



1/2 



(14) 



(15) 



The fit quality is demonstrated in Fig. [31 The root- 
mean-square relative error is 0.4%, the maximum error 
1.1% takes place at F = 10 and x = 1.67. It is seen that 
the fit is accurate for a classical Coulomb liquid (any 
F > 1) at least at a: < 3. This is sufficient to calculate 
the plasma screening enhancement of reaction rates (Sec. 
IIVB|) . It is remarkable that the fit is valid also in the 
crystalline solid (F = 400 and 800) although we have not 
included the crystalline data into the fitting. Clearly, 
H{r) in the solid and strongly coupled liquid (F 3> 1) is 
nearly the same, being "frozen" (almost independent of 
F) . This allows us to expect that Coulomb tunneling and 
the fusion reaction rate should not undergo significant 
changes when the temperature drops below the freezing 
temperature T„i. 



B. Enhancement factor 

Having H{r) we can introduce the mean-field reaction 
rate 



R 



MF 



pk 



7r/i(fcsr)3 



X / dE exp 

' Emin 



E 



P{E) 



, (16) 



where E is an energy of relative motion of colliding nuclei 
(with a minimum value i^min at the bottom of the po- 
tential well), exp(— iJ/fc^T) comes from the Maxwellian 
energy distribution of the nuclei, and S'pk is the S'-factor 



corresponding to the energy E at which the integrand 
has maximum. Finally, 



P{E) 



2V27Z 



dr 



H{r) - E (17) 



characterizes the penetrability of the Coulomb barrier at 
an energy E] rn and rj are classical turning points which 
are zeros of the expression under the square root; we can 
set r„ — > (an exact r„ should have been determined 
by nuclear interactions which we neglect in the Coulomb 
tunneling problem). Here we use a radial WKB approx- 
imation with the mean- field potential HJr). A similar 
approach was used by Ogata et al. [itI . [23| . The main 
difference is that Ogata et al. numerically solved a radial 
Schrodinger equation which should be equivalent to the 
WKB integration under the conditions of study. Notice 
that Ogata et al. employed less accurate mean-field po- 
tentials and obtained, therefore, less accurate results as 
discussed e.g., in [l^, [1^1 (also see Sec. IIIC|) . The ap- 
proach equivalent to our but with a less accurate mean- 
field potential was used by Itoh et al. 19] (Sec. Ill C|l . 

Putting H{r) = in Eq. (|17p we reproduce the well 
known result P{E) = 211 Z"^ / [hv) for a pure Coulomb 
barrier (with v = yj2E / Taking energy integral 
in Eq. p6)) by a saddle-point method we come to the 
thermonuclear reaction rate without plasma screening, 
R^^{H = 0} = i?th, given by Eq. (0. 

In the WKB mean-field approximation the enhance- 
ment factor of the nuclear reaction rate is 



MF 



R^''{H}/R^'{{)} 



MF r 



(18) 



The calculation of from Eq. (fT^ reduces to the eval- 
uation of P{E) and R^^{H) from Eqs. ^ and ^ for 
a given H{r). We have performed the integrations over r 
and E numerically (beyond the saddle-point approxima- 
tion) in all integrals. 

These results are naturally restricted by the WKB and 
mean-field approximations. They neglect fluctuative na- 
ture of plasma microfields; deviations from spherical sym- 
metry and dynamical evolution of these microfields dur- 
ing a tunneling event; corrections to H(r) due to quan- 
tum effects in ion motion; deviations from the first-order 
one-dimensional (radial) WKB approximation. In the 
thermonuclear regime with strong plasma screening all 
these effects are not expected to be strong. 



C. Analytic fit 



To facilitate applications of our results we have fitted 
the values of the enhancement factor ^ calculated 
from Eq. by an analytic expression. First of all, we 
notice that at 1 < F < 200 the function /io(F) can be 
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approximated as 

with Ai = 2.7822, A2 = 98.34, ^3 = V^- Ai/y/A^ = 
1.4515, Bi ^ -1.7476, B2 = 66.07, = 1-12, and 
B4 = 65. The accuracy of this fit is the same as the 
accuracy of the best fit (b) in Table HH At T > 1 
it gives /io*(r) « 1.0346 r and remains accurate at 
200 < r < 600. In addition, it reproduces the correct 
Debye-Hiickel asymptote h^^{r) = y^T^/^ at T <C 1. 
Note that the functional form of Eq. (fTH)) was suggested 
in Ref. [2^ to approximate the free energy of OCP. 

We have fitted the values of the enhancement factor, 
calculated on a dense grid of values 1 < T < 200 and 
< C < 8, as F^F = exp(/ifitp), 

h%{T, c) = h^„\T) + /if*(r, c) = C(r), (20) 

with 

r-r/(i + aC + /3C' + 7C')'/', (21) 

a = 0.022, P = 0.41 - 0.6/r, and 7 = 0.06 + 2.2/r. 
The maximum fit error of is « 30%; it occurs at 

C = 0.4 and F = 200, at which the enhancement factor 
itself is enormously large, F^f ^ 10^°. The fit accuracy 
is illustrated in Fig. [2] where the dashed lines (right ver- 
tical scales) show the ratio of fitted and calculated WKB 
mean- field values of gMF(O) [same as the ratios of fitted 
and calculated values of fj^/]. We have checked that 
an extension of ( from 8 to 50 at 1 < F < 200 does 
not change the initial fit accuracy. Moreover, if C < 50, 
the fit remains sufficiently accurate up to F ~ 600. For 
instance, at F = 400 the fit gives the values of F^^ 
which differ from the calculated values within a factor of 
2, and at F = 600 - within a factor of 10^ If F > 1 
s-i^d C Ij the main difference between F and F in 
Eq. (|2ip is determined by the term in the denominator 
containing /3 w 0.41 (while the terms containing a and 
7 are relatively small). Neglecting a and 7 for a mo- 
ment and treating /^^^ as a small correction, we obtain 
hi « -1.0346 X 0.41FCV3 ~ -0.141FC^ which is very 
close to -(5/32)FC^ ~ -0.156FC2 given by Eq. (A) of 
Table Hill and discussed in Sec. Ill CI 

To summarize, our fit gives very accurate values of 
F^f for all possible values of F and ( at which the WKB 
mean-field approximation is valid (see below) and, actu- 
ally, in much wider domain. 




c 



FIG. 4: (Color online) The screening function /if^(r,C)/r 
versus at six values of F. Dots with errorbars are Militzer- 
Pollock points [of group (ii) in Table HV] . solid lines are our 
mean-field WKB calculations, and short dash-dot-space lines 
are given by our fit. Dotted lines show the lowest-order result 
of Jancovici 16], while short-dash, long dash-dot-space, long- 
dash, and dot-dasli linespresent the results of Alastuey and 
Jancovici (AJ 78, Ref. Itoh et al. (IKM 90, Ref. fisl|), 

Ogata et al. (Oil 91, Ref. [H), and Ogata (Ogata 97, Ref. 
respectively. 
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V. PIMC, MEAN-FIELD WKB, AND OTHER 
RESULTS 

A. Overall analysis 

After calculating we have used Eq. pTjl and de- 
termined (?mf(0). The results are shown by solid lines 
in Fig. [5] for the same five values of rj which were taken 
by Militzer and Pollock This allows us to directly 
compare the PIMC and mean-field WKB approaches for 
all PIMC points [all groups (i)-(iii) of data points in Ta- 
ble llVj . Open symbols (right vertical scale) show ratios 
of the Militzer-PoUock to calculated mean-field results. 
The errorbars are those as reported Q in the PIMC sim- 
ulations of gMp(O) (with the corrections mentioned in 
Sec. imp . The overall agreement seems very satisfactory. 
Large differences take place in pycnonuclear points. Very 
strong differences (/mp (0)/(j'MF(0) ^ 0.04 and ~ 10"'' 
occur for 77 = 1 at F = 400 and 600, respectively. We 
analyze all these results below. 



B. Data of group (ii) 

After calculating the enhancement factor , we 
have presented it in the form ([7]) and determined . 
Using then Eqs. ^ and p9)) we have calculated — 
j^MF _ fi^(Yy This function is not very certain and has 
been a subject of debates (Sees. Ill CI and III Gp . Solid 
lines in Fig. 2] show our calculated values of /ii(F, C)/r 
versus C for six values of F = 1, 2, 5, 10, 40, and 100. 
Short dash-dot-space lines are given by our fit expression 
rSec UVCl) . 

Furthermore, taking the contact probabilities calcu- 
lated by Militzer and Pollock in the points which belong 
to group (ii) (Table [IV} and using Eq. pT|) we have de- 
termined the PIMC values of hi{T,() for the same six 
values of F as in Fig. 0] at several values of (. These 
data are plotted in Fig. 0] by dots, together with numer- 
ical errorbars of Militzer and Pollock One can ob- 
serve a remarkably good agreement between the mean- 
field WKB and PIMC results for the data of group (ii). 
Strongest disagreement occurs at the lowest F = 1, where 
the function hi introduces small contribution into the 
plasma screening enhancement of the nuclear reaction 
rates. The existence of real disagreement at F ~ 1 be- 
tween the PIMC and the mean-field WKB results could 
be checked in future more extensive PIMC runs. If real, 
this disagreement could be attributed to a fluctuative na- 
ture of the plasma potential at F < 1 (where Coulomb 
coupling is not too strong and can allow noticeable fiuc- 
tuations of the plasma potential from its mean-field val- 
ues). This would indicate that the PIMC results are more 
accurate at F < 1 than the WKB mean-field results. 

We have also compared the Militzer-PoUock and our 
results with some other calculations of /ii(F, C) discussed 
in Sec. Ill Cl In particular, the dotted lines in Fig.[5]show 
the basic lowest-order expression (A) from Table IIIII ob- 



r=6.i, c=o.22, r=io''K 




100 200 300 400 500 

r [fm] 



FIG. 5: Effective mean-field Coulomb potential U{r) for re- 
acting carbon nuclei at p = 5 x 10^ g cm~'^ and five tempera- 
tures (logj^Q r[K] = 9, 8.5, 8, 7.5, and 7). Shaded strips show 
Gamow-peak energy ranges; dotted lines are Gamow-peak en- 
ergies; dashed lines are thermal energies fcsT measured from 
the bottom of U{r). 
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FIG. 6; (Color online) Temperature dependence of the car- 
bon burning rate R (left vertical scale) and the inverse carbon 
burning time — R/rii (right vertical scale) at p = 5 x 10^ 
g cm"'' for the temperature range T > Tp of thermonu- 
clear burning with strong screening and for lower T (down 
to T ~ O.lTp). The upper line is obtained with the screening 
enhancement function h = /io(r) given by Eq. (b) in Table 
im Other lines employ the total function h = ho + hi cal- 
culated by Alastuey and Jancovici (AJ 78, Ref. [2J|), Itoh et 
al. (IKM 90, Ref. [li]), Ogata et al. (Oil 91, Ref. [ll), and 
Ogata (Ogata 97, Ref. The solid line plots our WKB 

mean-field calculation. See text for details. 

tained by Jancovici [l^; the short-dashed, long-dashed, 
and dot-dashed lines display the expressions (B), (C), 
and (D) calculated, respectively, by Alastuey and Jan- 
covici [ISIj Ogata et al. [13], and Ogata These 
results are shown for those ranges of F and C for which 
they were obtained in the cited publications (Table HIT)) . 
In addition, long-dash-dot-space lines present the func- 
tion hi{T,() which corresponds to the results of Itoh et 
al. [iB] (discussed in Sec. Ill C() . We see that the Militzer- 
PoUock and our results are in good agreement with ear- 
lier predictions of Jancovici. and especially of Alastuey 
and Jancovici, and Ogata [l^, but they are in worse 
agreement with the results of Ref. [l3|. Calculations of 
Ref. [l^ are also accurate; some whirls of correspond- 
ing curves at small C occur possibly because of simplified 
approximation of the mean potential in the cited publi- 
cation (see Sec. Ill Cp . 

As a byproduct of our calculations in the mean-field 
WKB approximation, we have determined the energy 
-Epk of the peak of the integrand function in Eq. ; it is 
the Gamow-peak energy modified by the plasma screen- 
ing effects. We have also estimated characteristic energy 
widths of the Gamow peak (at the half width of the peak 
maximum). In Fig. [5] we show the effective total radial 
mean-field Coulomb potential U{r) — Z'^e^ jr — H{r) for 
^^C ions reacting in pure carbon matter at p = 5 x 10^ 
g cm~'^. At this p, the ion-sphere radius is a=98 fm. Nat- 
urally, U (r) has a minimum at r « 2a due to Coulomb 



coupling. We plot U{r) for five temperatures, T = 10^, 
10^•^ 10^, 10'^•^ and lO"^ K. With decreasing T, the min- 
imum becomes more pronounced and finally "freezes" at 
F > 100. The shaded strips in Fig. [5] show the Gamow- 
peak energy ranges and the dotted lines show iJpk. The 
dashed lines present the thermal energy level /cbT mea- 
sured from the bottom of U{r). 

The first two upper panels in Fig. [S] refer to the ther- 
monuclear reaction regime with strong plasma screening 
(T > Tp). The next two panels are for a colder plasma 
(T = O.SeSTp and 0.114rp, respectively), while the lowest 
panel is for a very cold plasma (T = O.OSGlTp), certainly 
in the zero-temperature pycnonuclear regime. When the 
temperature decreases, the Gamow-peak energy range 
becomes thinner (note the difference of energy scales in 
different panels!) and shrinks to lower energies, together 
with iJpk- In the three upper panels the Gamow peak 
range is still at > [belonging to continuum states in 
a potential U{r)]. The energies from this range are much 
higher than fc^T supporting the statement that the main 
contribution into reaction rates at sufficiently high T 
comes from suprathermal ions. In these cases, the under- 
lying mean-field WKB approximation can be adequate. 
In the forth panel, the lowest energies of the Gamow-peak 
range become negative (drop to bound states) although 
-Epk is still positive and higher than ksT. The mean-field 
WKB approach may be qualitatively correct but quan- 
titatively inaccurate. At the bottom panel the Gamow- 
peak energy range fully shrinks to bound-state energies 
and the formal Gamow-peak energy becomes lower than 
ksT. It is clear that the mean-field WKB approximation 
breaks down at these low temperatures, and the formally 
calculated i?pk is inaccurate. Therefore, it is natural that 
the mean-field WKB results diverge from the PIMC ones 
at low temperatures. This divergence is seen in Fig. [2] 
(at highest values of F, especially at F = 400 and 600 for 

77=1). 

As follows from the above consideration, the mean- 
field WKB approximation can be valid for ( < 1.6 — 3 
[T > (0.1 - 0.2)Tp]. This is further illustrated in Fig. [5] 
which shows the carbon burning rate versus temperature 
at the same density p = 5 x 10^ g cm~'^ as in Fig. [5l The 
upper line is obtained with the simplified enhancement 
factor Fscr = exp[/io(r)], where ho{T) is given by Eq. (b) 
in Table [Til It is seen to be a good approximation in the 
thermonuclear regime (T > Tp) but gives qualitatively 
wrong results just after T decreases below Tp. It is easy 
to check that in the thermonuclear regime with strong 
screening hi is indeed a small correction to Hq. 

However, adding hi and using Fgcr — exp(/io -|- hi) 
greatly helps extending the strong-screening thermonu- 
clear results to lower temperatures, down to T (0.1 — 
0.2)Tp. Below Tp the function hi is no longer small but 
becomes comparable to ho and crucial to get physically 
reasonable results. All other lines in Fig.[6]are plotted by 
adding hi calculated in the various approximations. The 
solid line shows our mean-field WKB calculations which 
are nearly identical to our fit and to the Militzer-PoUock 
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FIG. 7: Pycnonuclear temperature-independent carbon burn- 
ing rate versus the density parameter rs. Open dots show 
Militzer-PoUock points. Shaded strip presents theoretical un- 
certainties of other calculations and solid line is the optimal 
model among these calculations • Filled circles with error- 
bars in the inset give the ratio of the Militzer-PoUock rates to 
the predictions of the optimal (Salpeter and Van Horn, SVH) 
model. 



PIMC results in the displayed temperature range. The 
short-dashed line corresponds to h\ and obtained by 
Jancovici [l^ and Alastuey and Jancovici [24j} for C < 1-6 
[Eqs. (b) and (B) in Tables E] and |TTT] . The dash-dot- 
space line is plotted using the fit formula of Itoh et al. 
[19] (C < 5.4). The dotted fine in Fig. [HI shows the results 
of Ogata et al. [13,111 (C < 2, Tables HI] and Hm and the 
long-dashed line shows the PIMC calculations of Ogata 
[l8^ (also C ^ 2, the same tables). We see that adding hi 
makes the reaction rate at T ^ O.lTp almost temperature 
independent, as it should be in the pycnonuclear regime. 

It is remarkable that all cited results (except for Refs. 
[ill,!!!) obtained at T > (0.1 - 0.2) Tp in different tech- 
niques and using various simplified assumptions (Sec. 
Ill C|) give actually almost one and the same reaction rate 
(almost the same curve in Fig. |6]) reproduced by the 
mean-field WKB approach. It gives us confidence that 
this approach is really valid at T > (0.1 — 0.2) Tp, and 
can now be considered as very reliable. It has been ex- 
pected by many authors (e.g., Rcf. [19]), and it is strictly 
confirmed by the Militzer-PoUock PIMC calculations ^ . 

C. Pycnonuclear Militzer-Pollock data [group (iii)] 

Finally let us analyze six Militzer-Pollock data points 
which correspond to temperature-independent pycnonu- 
clear burning [the data of group (iii) in Table IIV| . In 
Fig. [7| we plot the zero-tempcrature pycnonuclear car- 
bon burning rate as a function of the density parameter 
rs defined by Eq. ([3]). For simplicity, in this figure we 
use a constant (energy independent) astrophysical factor 
S{E = 1 MeV) = 3.2 x 10^^ barn MeV'^ for calculat- 
ing the reaction rates. The six Militzer-Pollock reaction 
rates are shown by open dots. 



We have compared these data with theoretical calcu- 
lations of other authors. The shaded strip shows other 
theoretical predictions taking into account the uncertain- 
ties of various approximations. It is plotted using the ex- 
pressions of Casques et al. ^10,] who analyzed calculations 
of different authors. The strip is restricted by the min- 
imum and maximum allowable reaction rates suggested 
in [i3] . The solid line is the optimal theoretical reaction 
rate, which is the static-lattice model of Salpeter and Van 
Horn 8] for a bcc Coulomb crystal. 

The consistency of the Militzer-Pollock data with other 
results is satisfactory. One should take in mind that the 
PIMC calculations in the pycnonuclear regime can be 
not too accurate because they require best computer re- 
sources. Moreover, the theoretical predictions of Ref. [T^ 
refer to pycnonuclear reactions in Coulomb crystal while 
the actual state of ions in the PIMC runs is unknown 
(not reported in 9 ]). As pointed out by many authors 
(see, e.g., Refs. [32.'33l[3^). strong zero-point vibrations 
of ions at > 90 — 160 prevent their crystallization even 
at T = 0. Therefore, three low-r^ points of Militzer and 
Pollock in Fig.[7]should correspond to pycnonuclear burn- 
ing in cold quantum liquid. However they do not deviate 
strongly from the predictions for Coulomb crystals. 

Nevertheless, the agreement of the pycnonuclear 
Militzer-Pollock points with the best theoretical predic- 
tion of Salpeter and Van Horn Q is not perfect. Filled 
dots in the inset in Fig. [7] display the ratio of the Militzer- 
Pollock to the best Salpeter- Van Horn reaction rates; the 
errorbars are those reported in (with the corrections 
mentioned in Sec. IIIip . We see that the Militzer-Pollock 
rates are noticeably lower. The nature of this difference 
is unknown. 

We expect that the Militzer-Pollock calculations in the 
pycnonuclear regime are not superior over other theoret- 
ical predictions at low T. More extended PIMC studies 
in the pycnonuclear regime would be helpful to reduce 
current theoretical uncertainties of the reaction rates. 

VI. CONCLUSIONS 

We have analyzed the recent Path Integral Monte 
Carlo (PIMC) calculations by Militzer and Pollock ^ of 
contact probabilities of atomic nuclei participating in fu- 
sion reactions in dense matter. We have compared these 
calculations with other theoretical predictions. In partic- 
ular, we have used a simple model based on WKB radial 
Coulomb tunneling of the reacting nuclei in the static 
mean-field plasma potential created by plasma ions. We 
have employed accurate Monte Carlo (MC) calculations 
of the mean-field plasma potential for a one-component 
strongly coupled plasma of ions and proposed a simple 
and accurate analytic fit to the plasma potential (Sec. 

Our main conclusions are as follows: 

1. We have found a very good agreement of the 
Militzer-Pollock PIMC results with the mean-field 



13 



WKB calculations for T > (0.1 - 0.2)Tp, i.e., in 
the thermonuclear regime, intermediate thermo- 
pycnonuclear regime and at highest temperatures 
of pycnonuclear burning (Table m. These results 
show good agreement with theoretical predictions 
of many authors (e.g., [H, E d 111] ) and can be 
considered as well established. 

2. There is a tentative slight disagreement of the 
PIMC and mean-field WKB results in the case of 
moderately strong ion coupling F < 1 but it cannot 
strongly affect the reaction rates. 

3. We have obtained a very accurate fit (Sec. IIV C|) to 
the plasma screening enhancement factors calcu- 
lated in the mean-field WKB approximation. The 
fit reliably describes the PIMC results in a wide 
temperature range T > (0.1 — 0.2)rp. 

4. New studies of Coulomb tunneling problem in the 
pycnonuclear regime are needed to obtain accurate 
reaction rates in this regime. 

The validity of the mean-field WKB method at F > 1 
and T > Tp could be expected (strong Coulomb cou- 
pling arranges quasi-order which may be well described 
by a radial mean- field without fluctuations). In contrast. 



there is little doubt that at T <C Tp the radial mean-field 
WKB picture is not true. At these low temperatures the 
reacting ions occupy quantum energy levels in their po- 
tential wells; they fuse along selected (anisotropic) close- 
approach trajectories [8] which is definitely beyond the 
mean-field radial WKB method. Because the accuracy of 
PIMC calculations 9] decreases at low T, a nice agree- 
ment between the PIMC ^9,] and mean-field WKB ap- 
proaches at T ~ O.lTp, which we formally reached, may 
indicate insufficient accuracy of the low-T PIMC results. 
New PIMC simulations would be most desirable to clarify 
this point. 



Acknowledgments 

We are grateful to B. Mihtzer, V. K. Nikuhn, and A. Y. 
Potekhin for useful remarks. Work of AIC and DGY 
was partly supported by the Russian Foundation for Ba- 
sic Research (grants 05-02-16245, 05-02-22003), by the 
Federal Agency for Science and Innovations (grant NSh 
9879.2006.2), and by the Dynasty Foundation. Work of 
HED was performed under the auspices of the US De- 
partment of Energy by the Lawrence Livermore National 
Laboratory under contract number W-7405-ENG-48. 



D. D. Clayton, Principles of Stellar Evolution and Nucle- 
osynthesis (University of Chicago Press, Chicago, 1983). [18 
P. Hoflich, Nucl. Phys. A 777, 579 (2006). [19^ 
T. Strohmayer and L. Bildsten, in Compact Stellar X- 
Ray Sources, edited by W. H. G. Lewin, M. Van der Klis [20 
(Cambridge University Press, Cambridge, 2006), p. 113. 

A. Cumming, J. Macbeth, J. J. M. in 't Zand, and D. [21 
Page, Astrophys. J. 646, 429 (2006). 

S. Gupta, E. F. Brown, H. Schatz, P. Moeller, and K.-L. [22 
Kratz, Astrophys. J. 662, 1188 (2007) 

D. Page, U. Geppert, and F. Weber, Nucl. Phys. A777, [23 
497 (2006). 

K. P. Levenfish and P. Haensel, Astrophys. Space Sci. [24 
308, 457 (2007). 

E. E. Salpeter and H. M. Van Horn, Astrophys. J. 155, [25 
183 (1969). [26 

B. Militzer, E. L. Pollock, Phys. Rev. B 71, 134303 
(2005). [27 
L. R. Casques, A. V. Afanasjev, E. F. Aguilera, M. 
Beard, L. C. Chamon, P. Ring, M. Wiescher, and D. G. [28 
Yakovlev, Phys. Rev. C 72, 025806 (2005). [29^ 

D. G. Yakovlev, L. R. Casques, A. V. Afanasjev, M. [30 
Beard, and M. Wiescher, Phys. Rev. C 74, 035803 (2006). 

E. E. Salpeter, Aust . J. Phys. 7, 373 (1954). [31 
D. G. Yakovlev and D. A. Shalybkov, Soviet Sci. Rev. [32 
Sec. E 7, 313 (1989). 

H. E. DeWitt, H. C. Graboske, and M. S. Cooper, As- [33 
trophys. J. 181, 439 (1973). [34 
Y. Rosenfeld, Phys. Rev. E 53, 2000 (1996). 
B. Jancovici J. Stat. Phys. 17, 357 (1977). 
S. Ogata, H. lyetomi, and S. Ichimaru, Astrophys. J. 372, 



259 (1991). 

S. Ogata, Astrophys. J. 481, 883 (1997). 

N. Itoh, F. Kuwashima, and H. Munakata, Astrophys. J. 

362, 620 (1990). 

A. Y. Potekhin and G. Chabrier, Phys. Rev. E 62, 8554 
(2000). 

H. DeWitt and W. Slattery, Contrib. Plasma Phys. 43, 
279 (2003). 

H. DeWitt and W. Slattery, Contrib. Plasma Phys. 39, 
97 (1999). 

S. Ogata, S. Ichimaru, and H. M. Van Horn, Astrophys. 
J. 417, 265 (1993). 

A. Alastuey and B. Jancovici, Astrophys. J. 226, 1034 
(1978). 

H. Kitamura, Astrophys. J. 539, 888 (2000). 

H. Kitamura and S. Ichimaru, Astrophys. J. 438, 300 

(1995). 

C. L. Jiang, K. E. Rehm, B. B. Back, and R. V. F. 
Janssens, Phys. Rev. C 75 015803 (2007). 

S. Ichimaru, Rev. Mod. Phys. 65, 255 (1993). 
J.-M. Caillol and D. Gilles, J. Phys. A 36, 6243 (2003). 
N. Itoh, N. Tomizawa, S. Wanajo, and S. Nozawa, As- 
trophys. J. 586, 1436 (2003). 

B. Widom, J. Chem. Phys. 39, 2808 (1963). 

D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 
(1980). 

G. Chabrier, Astrophys. J. 414, 695 (1993). 

M. D. Jones and D. M. Ceperley, Phys. Rev. Lett. 76, 

4572 (1996). 



